Multi-phase-field analysis of short-range forces between diffuse interfaces 
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We characterize both analytically and numerically short-range forces between spatially diffuse 
interfaces in multi-phase-field models of polycrystalline materials. During late-stage solidification, 
crystal-melt interfaces may attract or repel each other depending on the degree of misorientation 
between impinging grains, temperature, composition, and stress. To characterize this interaction, 
we map the multi-phase-field equations for stationary interfaces to a multi-dimensional classical 
mechanical scattering problem. From the solution of this problem, we derive asymptotic forms for 
short-range forces between interfaces for distances larger than the interface thickness. The results 
show that forces are always attractive for traditional models where each phase-field represents the 
phase fraction of a given grain. Those predictions are validated by numerical computations of forces 
for all distances. Based on insights from the scattering problem, we propose a new multi-phase-field 
formulation that can describe both attractive and repulsive forces in real systems. This model is 
then used to investigate the influence of solute addition and a uniaxial stress perpendicular to the 
interface. Solute addition leads to bistability of different interfacial equilibrium states, with the 
temperature range of bistability increasing with strength of partitioning. Stress in turn, is shown 
to be equivalent to a temperature change through a standard Clausius-Clapeyron relation. The 
implications of those results for understanding grain boundary premelting are discussed. 

PACS numbers: 68.08.-p, 61.72.Mm, 64.10.-|-h 
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I. INTRODUCTION AND SUMMARY 

Multi-phase-field models provide a powerful method to 
simulate complex interfacial patterns in a wide range of 
applications including polyphase and/or polycrystalline 
solidification 1 , 2, 4, H, @| , grain growth [7|, iSj] , as well 
as domain structures and solid-state phase transforma- 
tions . The equilibrium and non-equilibrium properties 
of isolated interfaces in multi-phase-field models are by 
now well-understood. Well-developed procedures exist 
for selecting model parameters in order to match some ex- 
perimentally specified set of interfacial free-energies and 
mobilities [5|, l6|. In contrast, the interactions between 
interfaces has remained comparatively more poorly char- 
acterized. Interfaces in phase-field models are inherently 
spatially diffuse. Hence, they interact when their dis- 
tance becomes roughly comparable to the interface width 
~ Those interactions can strongly influence the behav- 
ior of polycrystalline materials in many processes (such 
as sintering and solidification) where interfaces come into 
close contact at various processing stages. 

There have been a few studies of grain coalescence us- 
ing multi-phase-field models [l^, [IH as well as frame- 
invariant phase-field models with an order parameter rep- 
resenting the local crystal orientation [ij, [l^. Those 



studies have yielded useful insights but have been mostly 
numerical due to the inherent difficulty to treat analyti- 
cally the interaction between diffuse interfaces. 

In this paper, we develop an analytical approach to 
compute short-range interactions between diffuse inter- 
faces in multi-phase-field models. This approach is based 
on recasting the multi-phase-field equations for station- 
ary interfaces in the form of a classical mechanical scat- 
tering problem. A one-dimensional mechanical analog is 
standard for treating the properties of isolated station- 
ary phase-field interfaces It has also been used to 
treat interactions between nonlinear fronts in the real 
Ginzburg-Landau equation [T5|, which is analogous to 
the equation for a single phase-field. 

In a multi-phase- field context, the mechanical analog 
becomes higher dimensional, and hence more difhcult to 
analyze. It describes the motion of a point particle mov- 
ing in N dimensional space where N is the number of 
phase fields. In standard multi-phase-field models, each 
phase-field (jji describes the fraction of a given phase or 
grain orientation, which varies smoothly between zero 
and unity, with the physical constraint that X^i^i "Z*! ~ 1- 
The mechanical problem is therefore subject to this con- 
straint, but can also be formulated in iV — I dimensions 
after elimination of this constraint. This problem can 
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FIG. 1: Schematic plots of (a) disjoining potential Vex and 
(b) liquid film width W versus temperature for three qualita- 
tively different behaviors in an elementary material (i) purely 
repulsive, (ii) repulsive-attractive, and (iii) purely attractive. 
The solid (dashed) lines in (b) denote stable (unstable) equi- 
librium states. The disjoining potential Vex represents the 
excess interfacial free-energy due to the interaction between 
interfaces (i.e., the total excess interfacial free-energy minus 
twice the crystal-melt free-energy) and —dVex{W)/dW is the 
thermodynamic driving force causing interfaces to attract or 
repel each other. A uniaxial tensile (compressive) stress is 
predicted to have a similar effect as a temperature increase 
(decrease) in (b) through a Clausius-Clapeyron relation. The 
equilibrium state is monostable with a unique W below the 
melting point in (i) but can become bistable with two different 
W values with sufficient solute addition. 



be solved using conservation of total mechanical energy, 
which is the sum of kinetic and potential parts; the ki- 
netic energy is related to square-gradient terms in the 
multi-phase-field free-energy functional and the potential 
energy is just the bulk free-energy density in this func- 
tional, albeit with the opposite sign. The solution yields 
asymptotically exact analytical expressions for the forces 
between interfaces for distances W large compared to the 
interface thickness {W ^ (,)■ 

As a concrete example of application, we use the me- 
chanical analog supplemented by numerics to character- 
ize the interaction of crystal-melt interfaces. This in- 
teraction is relevant for understanding grain coalescence 
and grain-boundary premelting phenomena. The latter 
has been extensively studied experimentally [l6l [TtI ITsj 
and theoretically using lattice models [l^, |23|) atom- 
istic molecular dynamics or Monte Carlo simulations 
[U, m, (with earlier references therein), multi-phase- 
field models [l^, [llj , frame-invariant phase-field models 
[H, [H, and phase-field crystal models [H Even 
though grain boundary premelting is not fully under- 
stood, there is emerging concensus that it originates fun- 
damentally from a repulsive interaction between crystal- 
melt interfaces for high-energy grain-boundaries, which 
is directly relevant for the present study. This repulsion 
gives rise to the formation of an intergranular liquid film 
with a width that diverges at the melting point. 



The repulsive force responsible for the "premelting" of 
high-energy boundaries was computed in recent molecu- 
lar dynamics of pure Ni plj | and a two-dimensional phase- 
field crystal study of hexagonal crystals [25] . This force 
was found to decay exponentially with increasing dis- 
tance between interfaces in qualitative agreement with 
the form traditionally assumed in sharp-interface theo- 
ries [iQ. 26]. In addition, for low-energy boundaries, the 
phase-field crystal study revealed that this force is at- 
tractive at large distance but repulsive at short distance, 
as also recently observed in a molecular dynamics sim- 
ulation study of different grain-boundary types [2^ . In 
this attractive-repulsive case, the force vanishes at some 
intermediate equilibrium liquid film thickness, which re- 
mains finite at the melting point. For two grains with 
the same crystal orientation, in turn, the force between 
crystal-melt interfaces is purely attractive, in agreement 
with the fact that such grains generally coalesce to form 
a single grain. These three qualitatively different be- 
haviors: purely repulsive (i) for high-energy boundaries, 
attractive-repulsive (ii) for low-energy boundaries, and 
purely attractive (iii) for grains of the same orientation 
are depicted schematically in Fig. [TJ 

From a practical standpoint, liquid films can lead to 
a significant reduction of the shear resistance of a poly- 
crystalline mush and have been invoked recently to ex- 
plain hot cracking of metallic alloys during late-stage so- 
lidification [13, ll^. Therefore the ability to reproduce 
the correct cross-over from attractive to repulsive behav- 
ior with increasing grain-boundary energy is essential for 
modeling this phenomenon, and constitutes a stringent 
test for continuum models of polycrytalline materials at 
high homologous temperature. Ideally, a multi-phase- 
field model should provide enough flexibility to repro- 
duce force-distance {—dVex{W)/dW versus W) curves 
with the characteristics of Fig. [l] which can be com- 
puted from molecular dynamics simulations [U [2^ . 

A main finding of the present paper is that the stan- 
dard multi-phase-field formulation 0, d, @ is unable to 
reproduce the purely repulsive behavior for high-energy 
boundaries, corresponding to case (i) in Fig. [1] although 
it reproduces well the other behavior (ii), as well as (iii) 
that is a limiting case of (ii) for vanishing misorientation. 
In this respect, the standard multi-phase-field approach 
is more limited than the frame-invariant phase-field for- 
mulation that is able to reproduce all three behaviors 
[Tl, [13] . Using the mechanical analog, we show that this 
limitation of multi-phase-field models stems from the fact 
that phase fields generally represent phase (grain) frac- 
tions locally in space, and hence are constrained in the 
interval < < I in this interpretation. 

This result appears to contradict the finding of a purely 
repulsive behavior for a high-energy boundary in the re- 
cent multi-phase-field study of Cu-Ag alloys by Mishin et 
al. However, these authors used a modified multi- 

phase-field formulation that allows some phase-fields to 
become negative {(f)i < 0) in the region where diffuse 
interfaces overlap. Therefore, their results do not con- 
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tradict our finding that multi-phase-field models do not 
model pure repulsion when formulated in the traditional 
way where the 0i's represent positive phase/grain frac- 
tions. Rather, when interpreted in the light of the present 
analysis, the study of Mishin et al. [11] shows that the 
multi-phase-field approach can be modified to reproduce 
all desired behaviors in Fig. [T]with a less stringent phys- 
ical interpretation of the phase fields. 

In the present paper, we develop a different multi- 
phase-field approach where the free-energy landscape is 
inspired from the solution of the mechanical analog prob- 
lem. This approach abandons completely the interpreta- 
tion of the phase-fields as phase fractions and uses a min- 
imum number of phase-fields, as in a previous study of 
polyphase solidification (3|. This number is the same as 
the number of different grain orientations (two here for 
a bicrystal) , which is also the number of phase fields in 
the standard multi-phase-field formulation after elimina- 
tion of one phase-field using the constraint 4>i = ^■ 
The present formulation has the advantage of allowing to 
model the different interaction regimes in Fig. [T]by vary- 
ing a parameter that controls the sign and magnitude of 
the interaction between interfaces at large separation in 
an analytically predictable way. 

This formulation is developed first for an elementary 
material and then extended to a dilute binary alloy to 
investigate analytically and numerically solute effects on 
interface interactions. Solute addition is found to lead 
to the possibility of a qualitatively different behavior. 
Above a threshold concentration, two equilibrium states 
with different widths can coexist at some teniperature 
below the melting point, as also found in Ref. This 
temperature corresponds to a classical Maxwell point and 
the equilibrium state with larger (smaller) width is ther- 
modynamically stable (metastable) above this tempera- 
ture and vice versa below. We show analytically that 
this "bistability" follows from the fact that solute addi- 
tion makes the long-distance interaction between inter- 
faces more repulsive. Furthermore, we show that this 
effect becomes more pronounced for stronger partition- 
ing of solute between solid and liquid. Interestingly, this 
type of bistability was not observed in a recent atomistic 
study of grain boundary premelting in Cu- Ag alloys [23 | . 
In this study, the same boundary showed an attractive- 
repulsive behavior of type (ii) in Fig. [T]for both pure Cu 
and with Ag enrichment. However, this does not exclude 
the possibility of bistability for boundaries that already 
show repulsive behaviors in a pure case. 

Finally, we investigate the effect of uniaxial stress per- 
pendicular to a grain boundary on its premelting behav- 
ior. The coupling of solid-liquid phase change and stress 
is introduced by treating the liquid as a shear-free solid 
following the approach of Slutsker et al. [l^. Stress is 
shown to be equivalent to a temperature change through 
a standard Clausius-Clapeyron relation for a physically 
plausible choice of coupling between phase field variables 
and elastic energy, i.e. a tensile (compressive) stress cor- 
responds to heating (cooling). This prediction should be 



testable by atomistic simulations and experimentally. 

The paper is organized as follows. In the next section, 
we briefly review a simple sharp-interface model [13, HE] 
that provides an intuitive picture of repulsive and attrac- 
tive interactions. We then develop the mechanical analog 
in section IIIII We consider first the coalescence of two 
grains of the same crystal orientation that can be rigor- 
ously treated with one phase field. We then extend the 
approach to the more complex case of a bicrystal with 
three phase fields. We first discuss qualitatively why a 
repulsive behavior is difficult to obtain by examining the 
particle trajectories of the mechanical problem inside the 
Gibbs phase triangle and discuss how a simple triple- 
well one-phase-field model can produce pure repulsion. 
The mechanical analog is applied in section IIVI to com- 
pute explicit asymptotic forms of interaction forces for 
the multi-phase-field model of Ref. . The results con- 
firm the qualitative analysis of particle trajectories inside 
the Gibbs triangle. Next, in section [Vj we present our 
two-phase-field model of a pure bicrystal, which is a gen- 
eralization of the triple-well one-phase-field model, and 
analyze both analytically and numerically its properties. 
Solute and stress effects are then treated in sections IVII 
and lVIIl respectively. We conclude with a few remarks in 
section VIII. Technical details are given in several appen- 
dices where one appendix discussed difference between 
double-well and double-obstacle potentials. 

II. SHARP-INTERFACE THEORY 

The simplest picture of interface interaction is based 
on comparing at the melting point the excess interfa- 
cial free-energy of a dry grain boundary, 7gb, and the 
excess corresponding to two well-separated solid-liquid 
interfaces, 2"^ si- If 7g6 > 2,^sU the system can in principle 
lower its free-energy by forming a liquid layer, and the 
interfaces from two grains should repel each other. In 
contrast, if 7gf, < 2^sU the interfaces should attract each 
other so that the grain boundary remains dry. 

This picture can be extended to predict the width W of 
this liquid layer as a function of temperature by writing 
the total excess interfacial free-energy in the form 

AFe, = I^A/(T) + Ve^W) + 27./, (1) 

where A/ = fi — fs is the difference between the bulk 
liquid (/;) and bulk solid (/«) free-energy density and the 
sum of the other two terms represents the total excess 
interfacial free-energy. Close to the melting temperature 
Tm, 

Af{T)^L{T^TM)/TM, (2) 

where L is the latent heat of melting per unit volume. In 
addition, the quantity Vex{W) is the excess due to the 
interaction between solid-liquid interfaces, which can be 
assumed to have the simple form (Tol [26j 

VeAW) = i"fgb - 27./) exp(-T^/A), (3) 
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which interpolates between the Hmits of a dry grain 
boundary for — > and two well-separated solid-liquid 
interfaces for W — > -|-oo. The length A sets the range 
of the exponentially decaying interaction. As in recent 
studies 0, HH, we refer to Vex{W) as the "disjoining 
potential" by analogy with the disjoining pressure of fluid 
physics, i.e. the derivative —dVex/dW is the disjoining 
force that pulls interfaces a part when 7gf, > 2"fsi- 

This form reproduces the purely repulsive and attrac- 
tive cases (i) and (iii) in Fig. [T] when jgb is larger 
and smaller than 2"fgi, respectively. However, it does 
not reproduce the intermediate behavior (ii) with short- 
distance repulsion and long-distance attraction predicted 
in recent phase-field crystal 5^ and atomistic :23] model- 
ing studies. This limitation can be attributed to the fact 
that Eq. ([3]) assumes sharp interfaces and does not de- 
scribe the short-range repulsion associated with the for- 
mation of dislocations 5], which is still present for low- 
energy boundaries. While both the multi-phase- field [llj 
and frame-invariant [H, [l^] phase-field models also do 
not describe dislocations explicitly, the spatially diffuse 
nature of interfaces in those models suffices to produce 
qualitatively a short-range repulsion on a scale ^ ~ A and 
hence the intermediate behavior (ii). 

The temperature dependence of the liquid layer width 
is obtained by minimizing the excess free-energy given by 
Eq. ID) with respect to W, with Vex{W) given by Eq. 

. This miminization predicts a logarithmic divergence 
of as r approaches T^/ from below for jgf, > 2^ si , con- 
sistent with the behavior (i) in Fig. [Ijb). For < 2^ si: 
it predicts that the grain boundary remains dry over a 
finite superheated temperature range. The dashed line 
(iii) in Fig. [T{b) corresponds in this case to "unstable" 
equilibrium states. If interfaces are pulled slightly to- 
gether away from their unstable equilibrium separation, 
they attract each other until they join in the metastable 
dry grain boundary state with zero width. In contrast, 
if they are moved slightly apart, they repel each other to 
form a layer width of infinite thickness. 



III. MECHANICAL ANALOG 

A. Two grains of the same crystal orientation 

We consider first the coalescence of two grains with the 
same crystal orientation. A single phase-field is sufficient 
to distinguish between solid and liquid since both grains 
are equivalent. As depicted by case (in) in Fig. [TJb), 
crystal-melt interfaces are expected to attract each other 
for all separations W , since 7gb = 0.. As just explained 
at the end of the last section, this attraction implies the 
existence of unstable equilibrium states for T > Tm- The 
mechanical analog can be used to prove the existence of 
those states, and hence to conclude that the interaction 
is attractive. The free-energy per unit area of interface 




FIG. 2: Mechanical analog for coalescence of two grains of 
the same crystal orientation. The phase-field profile (a) cor- 
respond to the "coordinate" of a point particle moving in 
the "potential" U{<j)) = —fb{<l>) shown in (b) with x (the co- 
ordinate normal to the interface) measuring "time" in this 
analogy; /t is the bulk free-energy density. The trajectory 
is shown for a stationary, albeit thermodynamically unsta- 
ble, interface profile for T > Tm where the liquid has a 
lower free-energy than the solid (i.e., /i,(0) < /b(l) and hence 
(7(0) > (7(1)). In this case, the particle rolls down the po- 
tential energy landscape and then up to the turning point B 
after which it rolls back down and up to point A. This analogy 
also shows that a stationary interface profile cannot exist for 
T < Tm because of the absence of turning point in this case: 
the particle rolls past the liquid peak and never returns. This 
is consistent with the fact that interfaces from two grains of 
the same orientation cannot repel each other. 



has the form 



dx 



a f d(j) 



2\dx 



(4) 



where /fc(0, T) is the bulk free-energy density correspond- 
ing to a standard double-well potential with minima of 
equal height at T = Tm- A convenient form is 

/fc(0, T) = + 9T{<i^)L{T - Tm)/Tm (5) 

where fdw{4>) — h(j)'^{l — 0)^ has minima at and 1 cor- 
responding to liquid and solid, respectively, and gT{<l>) 
is a monotonously increasing function of with vanish- 
ing first derivative at and 1, and with 5t(0) — and 
5t(1) = 1. 

The equation for planar equilibrium solutions 
{5F/5(I) = 0) is 



dx"^ 



dU 



(6) 



where we have defined U = —fb- This equation has the 
form of Newton's law for a one-dimensional particle of 
"mass" (T and "coordinate" (j) moving in a potential U = 
—fb, with X measuring "time". The Hamiltonian for this 
dynamical system is the total energy, which is conserved 
in time. It is the sum H = K + U of the kinetic energy 
K = a{d(j)/dx)'^ /2 and potential energy U. 
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The proof of the existence of stationary solutions for 
T > Tm follows immediately form this mechanical ana- 
log. To see this, consider the phase-field profile corre- 
sponding to an unstable equilibrium solution for T > Tm, 
which is illustrated in Fig. [H^a). This solution depicts 
a situation where the attractive force between the two 
grains due to the overlap of the diffuse interface is bal- 
anced by the overheating that favors the liquid phase. In 
this analogy where x is time, the phase-field profile cor- 
responds to the trajectory of a particle in the potential 
[7, which has the form of a double-well potential turned 
up-side-down {U = —fb) with the liquid at a higher me- 
chanical potential energy (corresponding to a lower free- 
energy density) . The particle leaves the equilibrium point 
A, corresponding to the left grain, rolls down and then 
up the potential to reach the turning point B with zero 
velocity, corresponding to zero slope {d(j)/dx = 0) in the 
physical phase-field profile, and then rolls back down and 
up to the same equilibrium point A, which now corre- 
sponds to the right grain. It is clear that this A-B-A 
trajectory must exists as long as there is a turning point, 
which is always true for T > Tm- 

This mechanical analog can also be used to understand 
the divergence of W as the melting point is approached 
from above. For this, we note that the turning point 
approaches the liquid-peak of the potential energy as T 
approaches Tm- Therefore the particle will spend in- 
creasingly more time close to this peak as T becomes 
closer to Tm- Therefore, this time, and hence W in the 
analogy where time is x, must diverge as T ^ Tm- 

While this picture of the divergence is only qualitative, 
a quantitative understanding for large W is obtained by 
analyzing the trajectory close to the turning point and 
using conservation of mechanical energy. We sketch here 
the procedure and the details are elaborated in section 
IIVI Conservation of energy implies that 

H=^(^^J + UW^-LiT-TM)/TM, (7) 

where we have used the fact that the particle has zero ki- 
netic energy in the solid corresponding to the stationary 
point A in Fig. El and thus that H = U{1) = -fbil)- 
Applying this conservation law at the turning point cor- 
responding to point B in Fig. EJ we obtain that 

h^l^LiT-TM)/TM, (8) 

where we have used the fact that the value of at this 
point is small when the two interfaces are well separated 
and that 5t(0) yields a negligible contribution. This 
must be so because the turning point is physically located 
mid-way between the two interfaces. Since the phase-field 
decays exponentially in space away from the solid-liquid 
interfaces on both sides of this point, we would expect 
that 0tp « Aexp{-W/X) with A ~ ^, where ^ = 
is the interface thickness. This relation together with Eq. 
[5]predicts a logarithmic divergence of W as T—Tm — > 0+. 
Values for A and A are easily obtained by matching the 



solutions of Eq. ^ in the inner region close to the turning 
point and the outer regions close to the interfaces, which 
are both known analytically in this simple example. 

This analysis yields an analytical expression for the liq- 
uid layer width as a function of temperature, from which 
one can also obtain the disjoining potential using Eq. ([1]). 
For the present example, this yields 

VeAW) = -67,iexp(-M^/A), {W » ^), (9) 

with A = C/V2- The prefactor 67^; is three times larger 
than predicted by the sharp-interface theory, i.e. Eq. ^ 
with = (but depends on the precise definition of 
W for the diffuse interfaces). This is not surprising since 
the attractive interaction for large interface separation is 
governed by properties of spatially diffuse interfaces. 

B. Multi-phase-field model of a bicrystal 

The standard way to describe a system consisting of a 
liquid and two grains of different crystal orientations with 
a multi-phase-field model is to use one order parameter 
for each grain, chosen arbitrarily here as (/>! and (f)2 for 
grains 1 and 2, respectively, and a third (^3) for the 
liquid. In addition, the constraint 

(1)1 + h + h^'^, (10) 

is imposed consistent with the interpretation that each 0^ 
represents the volume fraction of the i*'' phase. This in- 
terpretation also implies in principle that < 0^ < 1 for 
each phase field but those constraints are not imposed. 
The range of variation of the phase fields depends gen- 
erally on the details of the free-energy functional. Stan- 
dard multi-phase-field models typically guarantee that 
< 4>i < 1 for all i. The same is true for the polyphase 
solidification model of Ref. [&], which is adapted to a 
bicrystal in section IIVI In contrast, in the formulation 
of Ref. [m , the ipi 's can become negative. In this sub- 
section, we restrict our attention to using a mechanical 
analog to draw general qualitative conclusions about in- 
terface interactions in a broad class of models where all 
phase fields vary in the interval zero to unity. 

The multi-phase-field free-energy functional can be 
written in the general form 

J [fkim, {V0,}) -I- MiM)] , (11) 

where fk is the "kinetic part" of the free-energy den- 
sity that contains gradient terms and ft is the bulk free- 
energy density. The former vanishes inside bulk phases 
while the latter remains finite. 

The stationary equations, which describe both stable 
and unstable equilibria, are given by 

5F 

— - Ao = 0, for i = 1,2, and3, (12) 
ocpi 
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grain 1 liquid grain 2 




liquid grain 1 



FIG. 3: Schematic representation of (a) the phase-field pro- 
files for a wet bicrystal and (b) the corresponding scattering 
trajectory (red dashed line) inside the Gibbs phase triangle. 
A physically admissible phase-field profile corresponds to a 
scattering trajectory where a particle leaves grain 1 with zero 
velocity and bounces from the liquid corner to arrive at grain 
2 with again zero velocity. The liquid corner is approached 
arbitrarily close as T — > Tm and W oo. If T > Tm, the liq- 
uid corner is at a higher mechanical potential energy than the 
corners corresponding to grains 1 or 2, and can therefore suc- 
ceed to produce this large angle hard scattering. In contrast, 
if T < Tm, the liquid corner is at a lower potential energy and 
the particle will generally scatter at a smaller angle from the 
horizontal axis, thereby leaving the Gibbs triangle. The exis- 
tence (absence) of trajectories for T > Tm {T < Tm) implies 
that the interaction between interfaces is generally attractive 
at large ly in a multi-phase-field formulation where trajecto- 
ries lie inside the Gibbs triangle. 



where 

is a Lagrange multiplier to satisfy the constraint (fTO|) . It 
is also possible to formulate the stationary equations by 
using the constraint (|10p to eliminate one of the phase 
fields, chosen arbitrarily here as (^3, directly in Eq. (fTTj) . 
The stationary equations then have, at least formally, a 
simpler form without constraint 

5F 

_ = 0, for i = 1 and 2. (14) 

For the effectively one-dimensional bicrystal geometry 
shown in Fig. [3l the stationary phase-field equations ([M]) 
are coupled ordinary differential equations with the in- 
dependent variable x. These equations are mapped to a 
classical mechanical problem for the motion of a particle 



TABLE I; Correspondence between the multi-phase-field sta- 
tionary equations describing a one-dimensional interface pro- 
file (Fig. |3]) and the classical mechanical problem of particle 
motion in a conservative potential where "dot" denotes dif- 
ferentiation with respect to "time" x in this analogy. 



Free-energy density 
./ = fk + fb 


Lagrangian 

L = T — U 


Free-energy F 


Action S 


Position X 


Time t 


Phase field (fii 


Coordinate qi 


Stationary phase field equations 
Scjiiix) 


Stationary action 
=0 


Generalized momenta 
df 

d(l>i 


Generalized momenta 
dL 

P' = "5^ 
oqi 


Hamiltonian 

H = ^Pi<i>i - / 

i 


Hamiltonian 
H = '^Piqi - L 

i 


Conservation law 
H ^0 


Energy conservation 
H = 



in a conservative potential by introducing the generalized 
momenta Pi = df /d(f)i, where we write a "dot" to denote 
d/dx to emphasize the analogy to classical mechanics. 
From those momenta, we can construct the Hamiltonian 

3 

H = Y,p,^,-f, (15) 

which is conserved in time [H — 0) and where f = fk + fb 
is the total free-energy density. Energy conservation also 
holds if the mechanical problem is formulated with the 
constraint pO)) since the latter is holonomic, i.e. it only 
depends on the phase fields and not their gradients and 
is "time" -independent. Both formulations without and 
with constraints are shown to be completely equivalent 
in Appendix |^ and we use here the formulation without 
constraint as described in Table HI 

Let us now examine the particle trajectories in the 
bicrystal geometry of Fig. [Sja). As in the last subsec- 
tion, the nature of the interaction for large separation 
{W :s> £_) can be deduced from the existence of parti- 
cle trajectories that correspond to physically admissible 
interface profiles close to melting. The interaction is at- 
tractive (repulsive) if stationary interface profiles exist 
for T > Tm (T < Tm)- It is useful to represent the 
particle trajectories in the standard Gibbs phase trian- 
gle shown in Fig. [Sj^b). The perpendicular distance of 
a point inside the triangle to an edge of the triangle is 
proportional to the volume fraction of the phase labeled 
at the corner opposite to this edge. Together with the 
constraint pO|) . this assigns a set of phase- field values 
{4>i,(l)2,4>3) for each point inside the triangle. 

The three corners of the Gibbs triangle correspond to 
minima of bulk free-energy density and hence to maxima 
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(a) 



grain 1 



X, 



grain 2 



(b) 

6o = 0, T < Ta/ 

eifi"! liquid e'!'"^ 



Liquid 

(scattering 

center) 




FIG. 4: Two-phase-field model with tunable interaction. The 
interaction is changed (a) by varying the distance (jio of the 
liquid free-energy minimum from the axis passing through 
the two solid minima. The scattering becomes "softer" , with 
the interaction switching from attractive to repulsive, as 0o 
is decreased. In the case (po = (b), the particle trajec- 
tory becomes simply one-dimensional and hops over the liq- 
uid minimum without being scattered. This trajectory can 
clearly only exist for T < Tm since U — —ft, showing that 
the interaction between interfaces is repulsive in this case. 



of the conservative potential U — ^ft- Consequently, a 
particle trajectory that connects the two grains, shown as 
a dashed line in Fig. ^h), leaves the grain- 1 corner with 
zero velocity at a; = — oo and ends at the grain-2 corner 
with zero velocity at a; = +oo. As T approaches Tm, the 
particle must approach the liquid corner arbitrarily close 
and spend a long time near that corner, corresponding 
to a large liquid film width. 

The remaining question is whether such particle trajec- 
tories exist in a slightly undercooled and/or superheated 
temperature range. For the one-dimensional mechanical 
analog of Fig. [2l the answer was clear since the point of 
closest approach to the liquid was a turning point. The 
particle only turned back if the liquid was at a higher 
mechanical potential energy, which required T > Tm 
since U = —fb- In the present case, the point of clos- 
est approach to the liquid (dark filled circle Fig. \^ is 
not a simple turning point since the particle has a fi- 
nite velocity at this point. Instead, the liquid corner 
acts as a "scattering center". A rigorous answer to the 
above question therefore requires a local analysis close to 
the liquid corner region to solve the scattering problem 
that connects incoming and outgoing particle trajecto- 
ries corresponding to diffuse solid-liquid interfaces. This 
analysis, described in section HVl for a specific choice of a 
multi-phase-field model, shows that scattering trajecto- 
ries inside the Gibbs triangle only exist above the melting 
point, and hence that the interaction between interfaces 
is always attractive for large separation. 

This answer can be qualitatively understood from the 
structure of the scattering problem with the help of Fig. 
[31 If r > Tm, the liquid corner is at a higher mechan- 
ical potential energy than the corners corresponding to 
grains 1 or 2, and can therefore succeed to scatter the 



particle at a large angle back towards the grain-2 corner. 
In contrast, if T < Tm, the liquid corner is at a lower 
potential energy and the particle will generally scatter 
with a smaller angle, thereby leaving the Gibbs trian- 
gle. The existence (absence) of trajectories for T > Tm 
(T < Tm) implies that the interaction between interfaces 
is generally attractive at large in a multi-phase-field 
formulation where trajectories lie inside the Gibbs trian- 
gle. 

This qualitative picture suggests how to construct a 
multi-phase-ficld approach to reproduce both attractive 
and repulsive interactions by relaxing the constraint that 
the trajectories lie inside the Gibbs triangle. The idea, 
which abandons the interpretation of the phase fields as 
volume fractions, is to construct a free-energy landscape 
where the free-energy density minima corresponding to 
the grains and the liquid are arranged in such a way that 
the scattering angle of the particle from the liquid cor- 
ner can be tuned to change the sign of the interaction, 
which is attractive for "hard" back scattering but repul- 
sive for "soft" forward scattering from grain 1 to grain 
2. A model with two phase fields cf) and ip built on this 
idea is shown schematically in Fig. [?1 and presented in 
more detail in section |Vl This model makes it possible 
to continuously change the interaction from attractive to 
repulsive by reducing the distance (p^ of the liquid free- 
energy minima from the axis passing through the other 
two solid minima. Reducing this distance reduces the 
scattering angle that vanishes for 0o = 0. In this ex- 
treme case, the liquid minima lies along the same axis as 
the two solid minima. Therefore the particle trajectory 
becomes simply one-dimensional and hops over the liq- 
uid minimum without being scattered. Such a trajectory 
can clearly only exist for T < Tm since U = — fb. This 
rigorously proves that the interaction between interfaces 
is repulsive in this limit of the model. 



IV. ANALYSIS OF MULTI-PHASE-FIELD 
MODELS 

In this section, we use the mechanical analog to 
compute analytically the large- distance interaction be- 
tween interfaces in standard multi-phase-field formula- 
tions where the particle trajectories lie inside the Gibbs 
phase triangle, as shown in Fig. O As discussed in the 
last section, this requires an analysis of the trajectory 
near the liquid scattering center, which corresponds to 
the liquid region between the two grains. We illustrate 
here this computation for the specific choice of the model 
of Ref . , but also consider other multi-phase-field for- 
mulations at the end of this section. The model of Ref. 
has two advantages for the present analysis. First, in 
the simplest case of equal interfacial energies, the phase 
field profiles are known analytically. Second, isolated in- 
terfaces between two phases, referred hereafter as binary 
interfaces, run exactly along the edges of the Gibbs trian- 
gle. Therefore the interface along a given edge does not 
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contain a spurious admixture of the phase labeled at the 
corner opposite to this edge (cj)^ — everywhere along 
the interface between grain 1 and grain 2, etc). We sup- 
plement our analysis by exact numerical computations of 
the forces for arbitrary distances between interfaces. 



A. Free-energy functional 



The individual contributions to the free-energy are 



i=l 



where h is a measure for the barrier height between the 
bulk states with the dimension of an energy density. The 
gradient energy is 



(17) 



where the parameter a plays the role of the mass in the 
mechanical picture. For the phase field model, h and 
a specify the solid-liquid free-energy and the interface 
thickness, see below. To allow for unequal solid-liquid 
and grain boundary interfacial energies, we add a grain 
boundary energy term which raises the free-energy well 
between the solid phases. 



fgb = ha4>i4>2{24)i4>2 + 303 



^1) 



(18) 



where only the dimensionless number a influences the ra- 
tio jgb /jsi 7 and b raises the free-energy bump only in the 
center of the Gibbs triangle but not along its boundary. 
We also introduce a coupling term 



(19) 



with the melting temperature Tm, the latent heat L and 
a thermal coupling function 



9t = — T- 



15(1-03)(1 

5)". 



(20) 



Then the total free-energy density is / = fdw + fgb + 
fk + fc, which is symmetric under exchange of (pi and 02 
as they represent the same solid phase only in different 
orientations. Hence /{, = fdw + fgb + fc in the above 
notation. The stationary equations are given by Eqs. (I14|) 
after elimination of 03 using the constraint that the sum 
of the phase fields equals unity. 



B. Liquid film width 

We now present a method to analyze the interaction 
between diffuse interfaces analytically and compare the 
findings to the numerical results. 



In the vicinity of the liquid corner, we can linearize the 
phase field equations and obtain 



(701 = 2h(j)i, i = 1,2, 



(21) 



where we have eliminated the third field. We note that 
the equations for both fields naturally decouple and have 
the same coefficients. This property will be discussed 
below in a more general context and become more trans- 
parent there. The total energy becomes in quadratic ap- 
proximation 



(16) H = -2h{<Pl + (f>l + (l)i(l)2) + Cr{:Pl + :j)l + <Pi,j)2) + L 



T-T, 



M 



Tm 



(22) 

On the other hand, the energy can be obtained from the 
limit X —^ ±oo, where it is i? = 0. The general solution 
of the linearized equations of motion is 



0i — Cii exp{Xx) + Ci2 exp(— Aa;) 



(23) 



with A = (2/i/ct)^/^. Symmetry with respect to exchange 
of the sohd fields demands Cn = C22 and C12 = C21, 
and then we obtain 



H = -4/i(4CiiCi2 + Ci\ + CI2) + L 



T - Tm 
Tm 



0. (24) 



For T < Tm one of the coefficients has to become neg- 
ative, but this immediately implies that the phase field 
coordinates will become negative at some moment. Since 
this contradicts the fundamental assumption that all 
phase fields have to stay in the range to 1, this shows 
that a solution with very wide liquid layer cannot exist 
below the melting temperature. Since we know that a 
solution must exist which connects to the macroscopic 
equilibrium solution W = 00 for T — Tm, this (unsta- 
ble) branch of solution must be located above the melt- 
ing temperature. Notice that there all coefficients can 
be positive, and we therefore cannot exclude the exis- 
tence of solutions. This is also illustrated in Fig. [5] for 
6 = 2, which we obtained from the numerical solution of 
the stationary phase field equations that were solved by 
a shooting method. Details of the solution procedure are 
described in Appendix [Bl Here we use the expression for 
the solid-liquid free-energy density 



(25) 



and the grain boundary energy Q 




-/gb = 2V2an I p{l-p)y/l + ap{l-p)dp. (26) 
/o 

In particular we see that the interfaces asymptotically 
attract each other, irrespective of the value of the grain 
boundary energy. 

For short distances W, however, the grains can also 
have a repulsive interaction, which of course does not 
contradict the asymptotic prediction. At sufficiently low 
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FIG. 5: Numerically computed liquid film width versus di- 
mensionless overheating for 6 = 2. The parts of the curves 
with negative slope are unstable for fixed temperature. The 
liquid layer thickness is defined as the distance between the 
point where the "solid" phase-fields cross the value 1/2. 



FIG. 7: Semi- logarithmic plot of the liquid film width as func- 
tion of dimensionless overheating. The numerically computed 
values (solid circles) are compared to the analytical prediction 
Eq. ((29]). We use 6 = 2 and -jgh/jsi = 2.2 here. 



grain 2 




0.001 



liquid 



grain 1 



FIG. 6: Numerically computed trajectories in the Gibbs 
phase triangle for different values of the overheating L{T — 
Tm) /hTu , which are labeled for each trajectory. The param- 
eters ^gb/'yst = 2.2 and 6 = 2 are the same as in Fig. [5] 



temperatures, the melt layer disappears, and the solution 
continues as a dry branch W — towards stronger under- 
coolings. The existence of these additional dry branches, 
which are not shown in Fig. [5] is a specific property of 
the model of Ref . [5] , and related to the absence of third 
phase contributions in a binary interface. 

For this model, the solution first runs nearby the edge 
of the Gibbs triangle which connects one solid phase with 
the liquid phase (see Fig. [6|). This is not mandatory for 
a general model, as many phase field models have third 
phase contributions in a binary interface, which implies 
that the trajectory deviates from the edge of the Gibbs 
triangle. We briefly discuss this case below. However, in 
the case of an interface between one grain and the liquid 
phase it is of course desirable not to have a contribu- 



tion of the other grain in the transition region, and much 
care was spent on fulfilling this requirement in the above 
model 

Since the binary solid-liquid interface profile is known 
analytically here, we can construct the solution for the 
asymptotic behavior for T — > T^: Approaching from 
X = — oo, let the solid- liquid interfaces be located at Xq = 
zLW/2 with a separation W/^ ^ 1, where 



1/2 



(27) 



is the thickness of an isolated interface with profile 



5l/2 



f 1 ± tanh ^ 
2 V 



(28) 



For ^ X ^ —Xq we can match its asymptotic behav- 
ior, 01 (x) ~ exp[— \/2(x — xo)/^] to the result from the 
linearization and obtain C12 — C21 — exp [-W/V2^]. On 
the other hand, the inner solution has to match asymp- 
totically a trajectory that passes along the edges of the 
Gibbs triangle for T — > . This implies that 0i /02 ^ 
for x — > 00 and therefore Cn — 0. Using the energy con- 
servation (23]) gives then 4/iCf2 = L{T-Tm)/Tm, which 
together with the above finding leads to the asymptotic 
behavior of the unstable branch 





■ L T 


— Tm 


\2h) 


Ah 


Tm 



for T ^ T 



M- 



(29) 



The comparison to the analytical prediction is shown in 
Fig. [71 Here we see explicitly that the liquid layer thick- 
ness diverges logarithmically when the melting point is 
approached, in agreement with the predictions of lattice 
models [l^ and molecular dynamics simulations [2l| . 
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FIG. 8: Numerically computed plots of disjoining potential 
versus liquid film width for b — 2 showing that the interaction 
is always attractive at large distance for the model of Ref. 01 • 
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FIG. 9: Semi-logarithmic plot of the excess potential as func- 
tion of the liquid layer thickness W . The data is compared to 
the analytical prediction Eq. (|33p . We use 6 = 2 here. 



C. Disjoining potential 

We can obtain the disjoining potential, using Eq. ([T]), 
which yields 



L{T ~ Tm) 



M 



(30) 



where /S.Fex{W) is the total excess free-energy, i.e. the 
total free-energy of the system minus the free-energy of 
a bulk solid phase occupying the same volume. The 
latter quantity is easily obtained by substituting the 
numerically computed phase-field profiles into the free- 
energy functional. The results plotted in Fig [51 con- 
firm the analytical prediction that the interaction is al- 
ways attractive for large W . The disjoining potential 
can also be predicted analytically by using the fact that 
d/^Ff,x{W) / dW = for a stationary interface, which 
yields the relation 



L{T - Tm) 
Tm 



(31) 



This relation reflects the fact that the grain attraction is 
compensated by the overheating. We therefore obtain by 
comparison with Eq. ([29]) 



W 



\2h) 



(32) 



which can be solved for VJ^ and integrated to yield 
VexiW) ~ -67,, exp (^-^^^ for W oo. (33) 

Here it becomes apparent that the grain boundary energy 
is not relevant for the long-range attraction. Instead, the 
strength of the interaction is solely set by 67s,, in contrast 



to the simple model As anticipated, the exponential 
decay takes place on the scale of the interface thickness 
^. This expression is compared to the numerical results 
in Fig. [9l showing an excellent agreement. 



D. Other multi-phase-field formulations 

Let us now briefly examine other multi-phase-field 
models than the one of Ref. ^5] with a more general ex- 
pression for the free-energy functional. We assume that 
the liquid phase-field is directly eliminated, so we do not 
need the Lagrange multiplier. Then the potential part of 
the free-energy is in quadratic approximation around the 
liquid point (f>i,4)2 = 0: 



fi 



b{(Pi, 



A/, 



and the kinetic energy is 



?l(fi2. 



(34) 



(35) 



Terms like 



cannot appear because they violate in- 



version symmetry. Positive definiteness requires a > |&|, 
and c > \d\, since the bulk liquid should be a stable so- 
lution; the exchange symmetry (pi ^2 is reflected by 
the above expressions. Notice that the above form of the 
kinetic energy contains also a case that is widely used in 
the literature (see e.g. Q and references therein) 



J2 ^ili^'^^^J 



(36) 



where the coefficients 7^ = jji are related to the inter- 
facial free-energies. In the vicinity of the liquid fixpoint 
they again reduce to terms of the above type. 
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The equations of motion are therefore 

c(f>i + d(l>2 — a(j)i+b<j)2, 
d(j)i+c(t>2 = b(j)i + a(j)2, 



(37) 
(38) 



We can define Ai = y/ {a + b)/ {c + d) and A2 = 
y/{a — b) /(c — d). The general solution is then 

(fii = ci exp(Aia;) + C2 exp(— Aix) 

+ C3 exp(A2x) + C4 exp(— A2X) (39) 

02 = ci exp(Aia;) + C2 exp(— Aia:) 

— C3 exp(A2a;) — C4 exp(— A2a;) (40) 

Symmetry requires C2 = ci and C4 = — C3, thus we have 

(pi — 2ci cosh(Aia;) + 2c3 sinh(A2a;), (41) 

(/)2 — 2ci cosh(Aia;) — 2c3 sinh(A2a;). (42) 



Then the Hamiltonian becomes 



H = -A/ - 4(a + b)cl + 4(a - b)cl = 0, 



(43) 



where the tilt term A/ corresponds as before to a devia- 
tion from the melting temperature (A/ < for T > Tm)- 
We can now distinguish three cases: Ai > A2, Ai < A2 
and Ai = A2. 

First, if Ai < A2, the "even" mode associated with cosh 
has the slowest decaying exponential. In the mechanical 
analog picture, it corresponds to a reflection of the parti- 
cle at the "liquid" potential hill. Notice that the match- 
ing constants Ci behave as q = q exp(— AiW^/2), with 
Ci being a number of order unity, which is determined 
from the matching of a single solid-liquid interface. This 
shows that in the first case C2 is exponentially small in 
comparison to ci in the limit W ^ 00 (r — > Tm)- Then 
this term does not appear in the energy conservation in 
this limit. Since a + b > the condition H = Q (which 
follows from the fact that the energy in the pure solid is 
zero) can only be fulfilled for negative A/, i.e. T > Tm, 
so the model is attractive at long distances. 

Second, for A2 < Ai, only the "odd" (repulsive) mode 
survives, so now ci is exponentially small compared to C2, 
and we can ignore the cosh part in the general solution. 
Then, however, the phase fields must become negative, 
which is forbidden. This mode corresponds to a particle 
that traverses the liquid bump, and therefore leaves the 
Gibbs phase triangle. 

Notice that in the case of unequal decay rates a pure 
binary interface cannot be free of third phase contribu- 
tions. If we assume that the solid with (pi — 1, (j)2 ~ 
for X — > —00 is in equilibrium with the melt ^3 = 1 
for X cxD, all growing exponentials must be suppressed, 
ci = C3 = 0, in the above general solution ([551 HO]) . Then, 
however, 02 — ±0i in the vicinity of the liquid fixpoint, 
which implies that a contribution of the other solid field 
(j)2 is always present. The equality of the exponentials is 
therefore a necessary condition for the absence of third 
phase contributions, as it is the case for the model above 



0; however, it is not a sufficient condition, and a coun- 
terexample is the model f3| , which is based on the kinetic 
energy expression ([551) . 

Finally, the model of Plapp and Folch 5] is a prototype 
of the last case of equal exponentials, Ai = A2. We re- 
frain here from performing a detailed general analysis of 
this case. Nonetheless we conclude that the mechanical 
analog, together with the restriction that the phase fields 
remain inside Gibbs triangle, poses a severe constraint. 
Therefore it is generally difficult to construct models that 
exhibit a long-range repulsion when the phase fields are 
interpreted as phase fractions. 

V. TUNABLE INTERACTION MODEL 

In this section, we present a simple two-phase-field 
model constructed around the idea that the large- 
distance interaction can be made repulsive by making 
the scattering trajectory of the particle softer in the me- 
chanical analog, as discussed at the end of section IIIII and 
illustrated in Fig. [H 

A. Model formulation 

A simple polynomial form of free-energy density with 
the structure of Fig. U] is given by 

h - h[{4^+l)'+^']-[{i:^l)' + (P^]-[^' + {cP+cPo)% (44) 

where 0o > measures the distance of the liquid mini- 
mum at (—00,0) from the axis passing through the two 
solid minima at ("0,0) — (~l:0)j (IjO)- A numerical 
example of the free-energy landscape is given in Fig. [10] 
and typical one-dimensional phase field profiles for a wet 
bicrystal are shown in Fig. [TTj In this model, the or- 
der parameter "0 has different values in the two grains 
of different crystal orientations {ip — ±1) and the liquid 
(-0 = 0), while only has different values in solid (0 = 0) 
and liquid (0 = — 0o). Therefore, it would be tempting 
to loosely interpret ^0 as a local measure of average crys- 
tal orientation, which vanishes in the liquid, and — 0/0o 
as a liquid fraction that varies from zero in the solid to 
unity in the liquid. However, such an interpretation has 
to be taken with caution for several reasons. Firstly, the 
model is not frame-invariant, hence the interpretation of 
■0 as a measure of local crystal orientation is not well- 
defined. Secondly, for 0o = 0, all the free-energy minima 
lie on the ip axis and [ipl can equally well represent a liq- 
uid fraction in this case. Thirdly, changing 0o has the 
same effect as changing the grain boundary energy and 
hence the misorientation. For these reasons, it is better 
to think as and "0 as the minimum set of two phe- 
nomenological order parameters necessary to construct a 
free-energy landscape with the desired properties. 
The gradient term (kinetic energy) is given by 



i(V0)2 + |(V0.)^ 



(45) 
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FIG. 10: Contour plot of the free-energy landscape for 0o — 
0.7, h = 1. The {ip = ±l,<j) = 0) minima correspond to the 
two different crystal orientations and the (■0 = 0,(j} — —(f'o) 
minimum to the liquid. For a — 1 the model is repulsive for 
00 < 0.69, in which case the scattering angle in the mechan- 
ical analog is sufficiently reduced for a particle trajectory to 
connect the two grains below the melting temperature. 



which has the desired property that, for finite 4>o, gi4>) 
varies from zero in the hquid to unity in the soHd and 
has vanishing derivative in the bulk in order not to shift 
the equilibrium values of the phase fields. Notice that 
treating the case where (po is exactly zero would in prin- 
ciple require a coupling function that also depends on 
-0. However such a complication is unecessary since we 
do not study this special case here, which was only dis- 
cussed in the context of Fig. |4] to motivate the model. 
The behavior of the model for 0o = and (po 1 are 
qualitatively very similar. 

The total free-energy density is then given hy f — fb + 
fk + fx, and the free-energy F is the volume integral of 
this expression. The stationary equations are 



0, 



5F 



= 0. 



(48) 



B. Liquid film width and disjoining potential 



0.5 



-0.5 



— 






t 

t 

» 





-4 -2 

X (h/a) 



1/2 



In the limit of an infinitely wide liquid layer, it is again 
sufficient to inspect the behavior in the vicinity of the 
liquid fixpoint ■0 = 0, = — 0Oj where the potential land- 
scape becomes to second order 



2\2„/,2 



+/i(l + 0^)^0^ (49) 
Then, the linearized equations of "motion" are 

cra0 = 1h{\ + (gf^, (50) 



= 2 



3 ^T~Tm 



Tm 



(0-f0o),(51) 



with the general solution 



FIG. 11: Phase field profiles for a one-dimensional wetted 
bicrystal geometry with a liquid layer sandwiched between 
two grains. Here a = 1 and 4>o = 0.7. 



where we introduce a > as additional parameter. In 
the mechanical analog, this corresponds to a tuning of the 
masses. Notice that we use the same parameters h and 
a as before; however, neither the interface width nor the 
interfacial free-energy can here be calculated explicitly, 
since the phase field profiles are not known analytically. 

Finally, to favor the liquid or solid states, we introduce 
a thermal tilt, 



fT = L 



T~Tj 



M 



Tm 



(46) 



corresponding to a homogeneous overheating or under- 
cooling with respect to the bulk melting temperature. 
We use the simple choice 



g(0) = l-(0/0o)'(3-f 20/0o), 



(47) 



ip{x) = c^+exp{\^x) + Cti,- exp{-\ti,x), (52) 
4>{x) = C0+ exp(A0a;) -I- C0_ exp{-X^x) - 00, (53) 



and growth rates 



2Ml + 0o) 



I2h{l 



6 j T-Tm 
4>l TKia 



(54) 
(55) 



Symmetry of the solution according to Fig. [TT] demands 
= —c^+ and C0_ = c^^. The Hamiltonian H — fj.— 
fb — fr becomes in quadratic approximation by energy 
conservation 



H = c^+CT [2aXl - 2A^(c0+/c.0+)^] = -L 



T-Tm 
Tm 



(56) 

Notice that in the case a = 1 the two exponential decays 
become equal when we approach the melting point, and 
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FIG. 12: Liquid film width as function of dimensionless over- 
lieating for a = 1.25. The parts of the curves with a nega- 
tive slope are unstable. In this model, the interface thickness 
^ = {(j/hY^^ and 7s! ~ {ahY^'^ so that the barrier height of 
the double well-potential h ~ jsi/£.- 



FIG. 13: Liquid film width as function of dimensionless over- 
heating for a = 0.75. The parts of the curves with a negative 
slope are thermodynamically unstable. 



this case will thus require some additional care. Let us 
therefore discuss a 1 first. 

For T Tj^j only the field with the slowest decaying 
exponential contributes, which is ip for a > 1 and (f> for 
a < 1. If the two sohd-hquid interfaces are far away from 
each other, they look (almost) the same as two single 
solid-fiquid interfaces, located at ±W/2, where W is the 
liquid layer thickness. We have asympotically e.g. for 
-W/2 < X < 0: ip ~ -*exp(-A^(a; -I- W/2)) and (f> ~ 
—(j)o + ^exp[—X^{x + W/2)], where both coefficients $ 
and ^E* are of order unity. Matching this to the above 
general solution ([5^ [55)1 gives then 

c^+ = -c^_ = *cxp(-A^VF/2), (57) 
C0+=c^- = $exp(-A0M^/2). (58) 

In the limit W ^ oo the weight factor c in front of 
the faster decaying mode (larger A) is exponentially sup- 
pressed in comparison to the other, and we can therefore 
drop its contribution in the Hamiltonian ((56)) . Hence, 
we can immediately conclude that the model is asymp- 
totically repulsive for a > 1, because solutions can ex- 
ist only for T < Tm, and vice versa for a < 1. No- 
tice that the value of (po is not relevant for this general 
long range interaction character; the asymptotic analysis 
makes predictions only for the limit T — > Tm (or equiv- 
alently — > oo), but the value (f)o can still significantly 
change the solutions with a liquid layer thickness of the 
order of the interface thickness. This behavior is shown 
in Fig. [12] for a — 1.25. Obviously, all cases are repulsive 
at large distances, but nevertheless a larger value of (po 
changes the pure repulsion and introduces a short scale 
attraction (with stable and unstable solutions above the 
melting temperature) and a first order transition char- 
acter. Here, we defined W as the distance between the 



points where tp crosses the values —1/2 and 1/2 respec- 
tively. 

The corresponding long-distance attraction is depicted 
in Fig. [T2]for a — 0.75, although the behavior is here less 
pronounced. 

To compare these results to analytical predictions, we 
need to determine 4* and $ (they are functions of (po 
and a) first; since, in contrast to the multi-order param- 
eter model of Ref. fE\ , the profile of a single solid-liquid 
interface is not known analytically, we have to find the 
matching constants numerically as follows: We set up a 
single solid-liquid interface at T = Tm, so the interface 
does not move. Assuming that the interface, i.e. the 
point ip = 1/2, is located at a; = and the solid phase in 
the domain a; < 0, we can look at the decay into the liq- 
uid region. For a; 3> 0, this decay is exponential, and we 
match it to ■)/) ~ 'I'exp(— A0x) and </) ~ $ exp(— A^a::)— (/)o, 
from which we get the desired prefactors \1/ and $. 

We can then extract the asymptotic behavior of the 
liquid layer thickness and the disjoining potential as be- 
fore. From the energy balance ()56p and the exponential 
prefactors ([571 [55)1 we get immediately 

r a^/^ja/hy/^ -L{T-Tm) . 

V2(l + 0g) ""A^^il + PlYTuh ""^"^ 

. 

The first asymptotic expression is of course applicable 
only for T < Tm, the second only above the melting 
point. Notice that here we had to evaluate (|55)) at 
T = Tm for the lowest order result. From this and the 
asymptotic relation V^^{W) ~ L{T — Tm)/Tm we get for 
the disjoining potential 

Ve^iW) ~ 2V2^^il + <l)l)a'/^ahy/^ X 
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FIG. 14: The ratio of the exponential prefactors as function 
of the liquid fixpoint position. For (po < 0.68 the model is 
repulsive, and attractive above this threshold. 



FIG. 15: Liquid film width as function of the dimensionless 
overheating, for a — 1. For (po < 0.68 the model is repulsive, 
and attractive above this threshold. 



X exp 



V2(l + (/)g) 
ai/2(CT//i)i/2 



W 



(60) 



for a > 1 and 



-2V2$2(l + 02)(^;j)l/2 X 

xexp\~-V2{l + (j)l)ia/h)-'^/^w] (61) 



for a < 1. 

Let us now look at the marginal case a = 1. There, at 
the melting point both exponentials have the same decay 
rate, and ip leads to repulsion, whereas (f) gives rise to 
attraction (which follows readily from the expression of 
the Hamiltonian), and it depends on the prefactors which 
effect is stronger. It will turn out, that the transition 
between attraction and repulsion is then controlled by 
(j)o, in agreement with the mechanical interpretation that 
the scattering angle for small (f>Q is small. 

We define the ratio of the exponential prefactors. 



2±t 



tp{x) 



exp[(A^ - X^)x] 



(62) 



Notice that this expression does not depend on tempera- 
ture in the limit T Tm, i.e. r = rQ + 0[{T — Tm)/Tm]- 
We can therefore determine the constant ro numerically 
from an isolated solid- liquid interface at T = Tm, and 
the result is shown in Fig. 1141 The reason why this 
is sufficient is that in the expression for the Hamilto- 
nian (|56p the common prefactor c^_^ is already of order 
(T — Tm)/Tm (as the right hand side), and therefore we 
need to evaluate the expression in brackets only at the 
melting temperature, i.e. to the order [{T — Tj^i)/Tj^j]" 
(exactly at the melting temperature the liquid layer is 
infinitely wide and therefore the exponential prefactors 
are zero). Then we get the solvability condition 

-T-Tm 



L- 



Tm 



Ah{l + ct>tYci^{l-rt) 



(63) 




W(h/0) ' 



FIG. 16: Disjoining potentials versus dimensionless film width 
for Q = 1 and three different values of 6n. 



Obviously, this equation has asymptotic solutions below 
the melting temperature only if Tq < 1, which is the case 
for |0o| < 0.68, and then the model is repulsive at large 
distances. The numerical results confirm this prediction, 
see in Fig. [TSl 

We can again calculate the asymptotic behavior ana- 
lytically, and obtain for the liquid layer thickness 



W 



[a/h) 



1/2 



In- 



-L{T-Tm) 



^/2(l + 02) 4(1 + 02)2vE,2(i _ ^l^TMh- 
Similarly, for the disjoining potential 



(64) 



K.(W^) ^2*%(l-r2)(a/i)i/2exp(-A^W^) (65) 

for W ^ oo. The full disjoining potential, as obtained 
from the numerical simulation, is shown in Fig. 1161 For 
the chosen parameters, the model is repulsive at short 
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FIG. 17: Semi-logarithmic plot of the liquid layer thickness as 
function of temperature (below the melting point). The data 
is compared to the analytical prediction Eq. (|64p . Parameters 
are a — 1 and (po — 0.5 here. 
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FIG. 18: Semi-logarithmic plot of the disjoining potential as 
function of the liquid layer thickness W . The data is com- 
pared to the analytical prediction Eq. (|65p . We use a = 1 
and 4>o = 0.5 here. 



distances even for a long range attraction (hard core re- 
pulsion). 

The analytical predictions are compared to the numer- 
ical results in Fig. [T7] for the liquid layer thickness and 
the disjoining potential in Fig. [15] for a = 1 and c/jq ~ 0.5 
(repulsive), which confirms the analysis. 

We therefore conclude, that the proposed model can 
describe both long-range attraction and repulsion. De- 
spite its simplicity the parameters can be tuned to cap- 
ture generic effects of many relevant materials. For 
a > 1, the model can also display bistability (coexis- 
tence) of "dry" and "wet" grain boundary states with 
different widths as shown in Fig. [T^lfor an intermediate 
value of 00- This bistability has also been predicted by 



a frame-invariant phase- field model of a bicrystal [12| . 
However, it has so far not been observed in molecular 
dynamics simulations of pure materials [2l|, [2^ . As 
shown in the next section, we find that solute addition 
can lead to bistability even for parameters of the model 
where bistability is absent in the pure limit. 



VI. SOLUTE EFFECTS 

We now extend the model to dilute alloys, correspond- 
ing to a phase diagram with straight solidus and liquidus 
lines. This dilute limit is described by adding to the free- 
energy density the contribution due to solute addition 



fc 



RT, 



M 



Wo 



(c In c - 



x(0)Aec, 



(66) 



where R is the gas constant, vq is the molar volume, and 
c is the mole fraction of solute assumed much smaller 
than unity. This contribution includes the standard en- 
tropy of mixing term and a partitioning term that distin- 
guishes between the energy density of impurities in solid 
and liquid via the coupling function x{4')- This function 
varies from in the liquid to 1 in the solid and may be 
chosen equal to g{4>)- A dependence on -0 can also be 
introduced to influence the segregation of impurities at 
the grain boundary, but we do not investigate this effect 
here. 

The concentration field obeys in equilibrium the con- 
dition 



5F_ 



II. 



(67) 



Since it enters the free-energy functional without gradi- 
ent terms (in the mechanical analog, the "coordinate" c 
belongs to a particle without mass, which follows the mo- 
tion of the phase fields "instantaneously"), we can elim- 
inate it and rewrite the phase field equations as derived 
from the grand potential. We obtain from the expression 
above 



where we defined 



exp 



exp 



RTm 



/ Vpfi \ 

\RTmJ 



(68) 



(69) 



Here, we immediately identify the meaning of the parti- 
tion coefficient fc, 



k — exp 



vaAe \ 
RTmJ ' 



(70) 



since we get for the concentrations of an (infinite) solid- 
liquid equilibrium system 



(71) 
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Notice that for a thin Hquid layer the concentration dif- 
fers there from the expression ([55]) . since the phase field 
(j) does not fully reach the liquid value —(po- 

We can change the ensemble and eliminate the con- 
served field and replace it by the intensive variable fi. The 
adequate thermodynamic functional is then the grand 
potential, from which the phase field equations can be 
derived variationally. This implies that the mechanical 
analog holds with the potential energy of the particle now 
determined by the grand potential density instead of the 
free-energy density. Then, the impurity contribution to 
the grand potential, lOc = fc^ l^^c, is 



RTm 



c(0). 



(72) 



Again, for the equilibrium of two bulk phases, the grand 
potential 



must be minimized, i.e. 



- '^^ - 

Sip ^ 5(j) ^ 



(73) 



(74) 



which implies that its density is equal in solid and liq- 
uid for W ^ Qo. Here, uj — uJq + fb + It + fk- Since 
/b = /fc = in both infinitely large bulk states, we get im- 
mediately Wc(0 = 0) + L(Teq - Tm)/Tm = Wc(0 = -0o), 
and therefore 



which describes the straight liquidus line with slope 



(1-fc). 



(75) 



(76) 



Expanding again up to second order around the liquid 
fixpoint we get 

= + l v-c|-^)(lnfc)(^+^o)^ 

TM{l-k) ' 2Tm(1-A:) ' AV^t^w; 

(77) 

with x" — x"{4> ~ ^4'^)- From the total grand potential 
we get the linearized equations of motion with — 0J~fk 




2h{l + (blr , 6 ^T,q - Tm 



(78) 
(79) 



mL 



Tm(1 - k)a 



c[^^'(lnfc)x" 



1/2 



(80) 
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FIG. 19: Liquid film width as function of overheating for dif- 
ferent solid concentrations. The addition of impurities makes 
the model more repulsive at large distances, and leads to a 
first order character at shorter liquid layer thicknesses. Pa- 
rameters are a = 1.0, 00 = 0.5, mL/TMh — —1, k — 0.5. 



The solution for the linearized phase fields has again the 
structure ([5^ [55]) . Obviously, the decay rate of the liquid 
field 4> is modified in comparison to the pure case, and it 
becomes larger here. This means that the model becomes 
more repulsive through alloying, and the effect is more 
pronounced for stronger partitioning. For the particular 
choice X — 9i with g being the thermal coupling function 
(|T7)) we obtain 



a Tm 



Tm 



Infc 
1 - k 



1/2 



(81) 

which has to be compared to the decay rate of ij) given by 
Eq. The influence of impurities is shown in Figs.fTOl 

and [20] as function of the temperature deviation from 
equilibrium and fixed chemical potential; for all numer- 
ical calculation x = g is used. For convenience, the 
latter quantity is expressed through the concentration in 
the solid far away from the interfaces. It is equivalent to 
the notion of the chemical potential through the relations 
and ([7T|) . In both cases, the addition of impurities 
leads to a pronounced first order character, and an en- 
hanced repulsion at large distances with higher impurity 
concentration. For the particular marginal choice a — 1 
in Fig. [201 which is attractive in the pure case by the 
choice of 00; the model becomes immediately repulsive 
through the presence of impurities. 

Fig. [21] shows the profiles of the phase fields and -ip 
and the concentration c as function of the position for pa- 
rameters as in Fig.[Tl = 0.25 and L(T --Teq)/TMh = 
—0.02. Here, two stable and one unstable solution exists, 
which differ by the melt layer thickness. For the solution 
with the widest liquid layer a rather pronounced liquid 
phase exists, i.e. the phase field (p is almost stationary in 
the center, but it does not fully reach its bulk equilibrium 
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FIG. 20: Liquid film widtli as a function of dimensionless 
overlieating. Wfiereas tlie pure case is attractive, the addition 
of impurities leads to a long-range repulsion. Parameters are 
a = 1.0, </)o = 0.8, mL/Tuh = -1, fc = 0.5. 
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value "(pQ. Consequently, also the impurity concentra- 
tion is significantly larger than in the bulk solid. For the 
solutions with the thinner width, the grain boundary is 
almost dry, and the concentration only slightly increased. 

Again, for x — 9i starting with an attractive situation 
with a < 1 without impurities, the long-range interaction 
becomes repulsive for 

(82) 

because then the decay lengths of the "attractive" field 
(j> and the "repulsive" field iJj become equal at the co- 
existence point; the solution of this equation defines a 
critical temperature T*(a,(j3o,m,k). We can then de- 
fine a (dimensionless) deviation from this value as t — 
{T* — T)/T*, and the numerical results are shown in 
Figs, mi and [231 Here, we keep the temperature con- 
stant and vary the chemical potential, and the behavior 
is qualitatively similar to the curves with fixed chemi- 
cal potential and varying temperature. The equilibrium 
chemical potential ^eq is given by the expressions (|69p . 
(ffTj) and ([75)) . Here, from the given temperature T the 
equilibrium chemical potential /ieq can be calculated; a 
change of the "supplied" composition in the solid phase 
far away from the grain boundary allows then to vary 
the chemical potential. The results confirm the analyti- 
cal prediction that the system becomes repulsive at long 
distances below the critical temperature T* . 

In the same way as before, we can calculate the asymp- 
totic energy balance using the fact that the Hamiltonian 
H ^ fk ~ Wp is constant, which yields 

- c[''\l - k)6fi - 2cl+aaXl ~ 2cl+aXl,, (83) 
with 5fi — ^ — ^eq- The logarithmic divergence of the 
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FIG. 21: Profiles of phase fields and impurity concentration 
for a temperature inside the bistable range of Fig. 1191 Cs — 
0.25 and L(T - Teq)/TMh = -0.02. For this temperature, 
three solutions exist with W decreasing from top to bottom. 
The top and bottom solutions can coexist whereas the middle 
solution is thermodynamically unstable. 



liquid layer thickness at /igg follows again immediately 
from the preceding relations, and we can take into ac- 
count also the effect of the second exponential, resulting 
in the implicit relation 

(1 ~ k) ^ 

-2a$2cxp(-A^^,iy)A2,,), (84) 

where the matching constants '^(a,(j)o,m,k,T) and 
$(«, 4>Q, TO, k, T) are determined from a single solid-liquid 
interface at bulk equilibrium as before. For large sepa- 
ration, of course only the slowest decaying exponential 
contributes, but the inclusion of the next term can lead 
to a substantial better agreement with the numerically 
obtained result, as shown in Fig. [24l 

Finally, we checked the influence of the partition co- 
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FIG. 22: Liquid film width as function of tlie cfiemical po- 
tential deviation from bulk equilibrium for different values of 
the dimensionless reduced temperature t (see text). The inset 
shows a magnification around the origin to demonstrate the 
transition from an attractive to a repulsive behavior if the 
temperature is lowered. Parameters are a — 0.9, (po ~ 0.7, 
mL/TMh = —1 and k = 0.5. 



6 




1 I • • • 1 

-0.2 -0.15 -0.1 -0.05 

{^-Heq)/h 



FIG. 23: Liquid film width as function of the chemical po- 
tential deviation from bulk equilibrium for different values 
of the dimensionless overheating. Parameters are a — 1.25, 
(l>o = 0.75, mL/TMh = -1 and k = 0.5. 



efficient on the results. Fig. [21] shows the melt layer 
thickness as function of temperature for fixed chemical 
potential. We see in general that a stronger partition- 
ing system can be overheated more. In particular, for 
the example shown here, no equilibrium solution exists 
above the melting point for k ~ 0.5, whereas alloys with 
smaller partition coefhcient exhibit stable and unstable 
solutions also above the melting point. 

We have only studied here a free-energy functional 
without a gradient term ~ (Vc)^. With the inclusion 
of such as term, the equilibrium condition (|67p has the 
same form as the phase field equation Therefore, 



FIG. 24: Comparison of the numerical data for the deviation 
of the chemical potential Sfi versus the melt layer thickness 
W with the asymptotic prediction for a binary alloy. The pa- 
rameters are a — 0.9, (po — 0.7, k = 0.5, mL/TMh = —1 and 
L{T — Tm) /TKih — —0.45. At large distances the interfaces 
repel each other, and the contribution of the second expo- 
nential leads to an attraction. Asymptotically, the behavior 
is determined by the slowest decaying exponential only; for 
smaller separations, the behavior is well described by the con- 
tribution from the two slowest decaying modes. 
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FIG. 25: Liquid film width as function of overheating for a 
repulsive case and different partition coefficients. The chem- 
ical potential is the same in all cases, expressed through the 
equilibrium liquid concentration c^'^''^ — cf'^'' /k — 0.5. The 
other parameters are a = 1, tpo = 0.5 and mL/TMh = — 1. 



the same type of analysis that exploits a mechanical ana- 
log can be performed with a concentration field that now 
possesses a "mass" . 

In summary, a pure system that is attractive at melt- 
ing, can become repulsive with solute addition above 
some threshold concentration. Furthermore, solute ad- 
dition can lead to bistability (i.e., existence of stable 
and metastable states at the same temperature on ei- 
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ther side of a Maxwell point) with the effect being more 
pronounced for stronger partitioning (smaller k). A nu- 
merical estimate of the temperature range of bistabil- 
ity is useful to examine if first order hysteretic transi- 
tions between liquid films of different widths could be 
observed. In Fig. [251 bistability extends over a dimen- 
sionless temperature range \L{T — T^q) /{Tmh)] of almost 
0.1 for k = 0.01. Since h ~ Isi/^,, we obtain that 
\T-T^q\/TM ^ 0.l7si/(LC)- With ^ - 1 nm and typical 
values of L and 7^; for metallic systems (e.g., L « 3 x 10^ 
J/m'^ and 7^; « 0.3 J/m^ for pure Ni), we obtain that 
\T—Teq\/TM ~ 10~^. So, according to the present model, 
bistability should be present over a temperature range 
below melting of the order of tens of degrees. 

VII. STRESS EFFECTS 

Stress effects can have a strong influence on microstruc- 
tural evolution in the presence of dry and wet grain 
boundaries at high homologous temperature, as in the 
practical case of hot cracking of metallic alloys [l3] . Here 
we limit our study to the effect of a uniaxial stress ap- 
plied along an axis perpendicular to the grain bound- 
ary. We couple phase change and elasticity by modeling 
the liquid as a solid with vanishing shear modulus, as 
in previous studies of solidification under stress HI], the 
Asaro-Tiller-Grinfeld instability [3]| . and fracture asso- 
ciated with a phase change [30|. In this approach, the 
local displacement vector (with two degrees of freedom in 
two dimensions) is also represented in the liquid, whereas 
the local state in this phase is fully characterized by the 
hydrostatic pressure. The remaining degree of freedom 
allows to represent the slip of a liquid at the solid-liquid 
interface, as an incoherent boundary between two solids. 
For the present purpose of a one-dimensional analysis, 
the situation is even simpler since we assume that the 
liquid film cannot expand in the direction parallel to the 
interface. Then solidification shrinkage due to a density 
difference between solid and liquid can be rigorously de- 
scribed as an "eigenstrain" . 

Since the focus is here on interactions on short scales, 
the issue of proper coupling of the elastic fields to the 
local phase field arises. As before for the temperature 
and impurity coupling, this relationship is not unique, 
and different choices can lead to the same sharp inter- 
face limit. Therefore, additional physical assumptions or 
input from other sources is required. We discuss here two 
different choices of the coupling function to illustrate the 
consequences of this effect. 

The first choice for the additional free-energy contri- 
bution is given by 

/ei = ^Ae|j -|-/2e|fe (85) 

with 

e., =ey-[l-.9('/')K^f- (86) 



The strain is derived from the displacements Ui as 

Here, we assign the eigenstrain e'^' to the liquid phase, 
and we use the same coupling function g{(f>) as above. 
Also, we define averaged elastic constants 

A = (7(0)A, + [l-5(0)]Ai, (88) 

where As and A; are the first Lame coefficients of the 
solid and liquid phase, respectively. Similar definitions 
are used for the shear modulus p,. The entire free-energy 
F depends then also additionally on the displacement 
field, and equilibrium requires that F is also minimized 
with respect to this new degree of freedom, which implies 
static elasticity, daij/dxj — 0. The stress tensor is here 
given by Hooke's law, 

aij = 2/2ey -I- Xlkk^ij- (89) 

The advantage of this model is that in the case of equal 
elastic constants in a one-dimensional case the influence 
of stress can be mapped to a temperature tilt: We assume 
/I = (since the liquid phase has no elastic response to 
shear and we demand the equality of the elastic constants 
in both phases), A^ = A; = A, and the only nonvanishing 
displacement component Ux depends only on x. Then 
elastic equilibrium, SF/Sui — 0, requires that the stress 
a^x is spatially constant, with 

Txx = Xsl^xx - (l - 9m4?]- (90) 
The elastic free-energy density becomes then 

fel - ^^L/A.. (91) 

Notice that the proper underlying boundary conditions 
for a minimization of the total free-energy are fixed dis- 
placements; although the stress is spatially constant it 
varies if the interfaces move, and we get the additional 
contribution to the equations of motion 

^=-..ff'(0)4°\ (92) 

which has exactly the same structure as the driving force 
term arising from the thermal tilt (|46p . We can then 
immediately identify the driving force term in the phase 
field equation due to elasticity with the one which stems 
from thermal effects if we relate 

T - Tm _ (0) 

= eli'axx, (93) 

which is the classical Clausius-Clapeyron relation since 
eii is the relative volume change , i.e. e^^ = Av/v. As a 
result, the application of a stress is equivalent to a change 
of temperature, and the abscissae of all previous plots 
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FIG. 26: Liquid film widtli as function of tlie imposed fixed 
displacement 5 at solid surfaces at melting (T — Tm)- Param- 
eters are a = 1, L = 100 {a/hy^^, e[°> = 0.1 and \/h = 500. 



of liquid film width versus dimensionlcss overheating in 
section IVl can be relabeled with {Av/v)axx/h in place 
of L{T — Tm)/ {Tuh). For an applied stress to change 
the liquid film width appreciably, it should produce an 
equivalent temperature change of at least a few degrees. 
For a typical relative volume change of a few percent and 
L a few times 10^ J/m^, this stress axx = L{Av/v)~^{T— 
Tm)/Tm should be in the range of hundreds of MPa, and 
therefore of the same order as a typical yield stress. 

Apart from situations with given stress, we can also 
consider the case of given displacement, and we can ob- 
tain it directly from the results for fixed stress. From 
Eq. ([90]) we obtain the total displacement 

L L 

S := J dx = + e|"^ j [I - (94) 



which depends now explicitly on the system size L. We 
note that here always a "macroscopic" stress free solution 
exists with txx = in the solid and txx = in the liquid 
(at T = Tm). This imphes 

We^^ = S (95) 

independent of the elastic parameters; this equation re- 
flects that the applied displacement is compensated by 
the volume change during melting, and the liquid layer 
thickness adjusts itself such that the elastic energy is min- 
imized (at T — Tm), i.e. the system is stress free. The 
above relation is of course a sharp interface prediction, 
therefore valid asymptotically for large S/e^\ 

The numerical results are shown in Fig. 1261 they 
exhibit the correct asymptotic behavior, which does not 
depend on the interaction of the interfaces, as it be- 
comes negligible for large interface separations. Notice 
that equilibrium liquid layer thicknesses for different sys- 
tem sizes look different if plotted versus the average strain 



5/L, but they collapse to the same curve if drawn as func- 
tion of stress, provided that the system size is much big- 
ger than the interface thickness, L 3> {a /Kf-I"^ . Also, it is 
worthwhile to mention that the stability of the branches 
has changed in comparison to a case with prescribed 
stress: Whereas for fixed stress the asymptotic branch 
with negative slope is unstable, solutions on the macro- 
scopic branch for fixed displacement are stable, since it 
corresponds to an energy minimum; nevertheless they 
correspond to the same solution. The same behavior will 
of course also occur for the pure thermal coupling if in- 
stead of the temperature the heat content is kept fixed. 

We note that the stress-temperature duality is a con- 
sequence of the chosen coupling function and the proper 
free-energy functional and not necessarily true for other 
choices, and it holds in general only for equal elastic con- 
stants. To make this effect more transparent, let us con- 
sider the special situation of a stress free state, where we 
get for the above coupling /gj = 0. An alternative way 
to implement the elastic energy is [2^ 

hi - \k + [1 - 9(.m^n ~ 4?)') ' (96) 

where we, for simplicity, directly assumed /i = and 
equal Lame coefficients. It differs from the above expres- 
sion only by the averaging in the interface region and has 
the same sharp interface behavior. Now we get from the 
mechanical equilibrium condition, 

Oxx = A, (e..5((/)) + (1 - .9(0))(e.. - ef )) = 0, (97) 

where we assumed again a stress free situation. Multipli- 
cation with (txx — ^ix) integration over x gives 

L 

= - J txx 9{^) dx (98) 



with Fei = J feidx, which does not vanish in a stress 
free situation. Notice that according to Eq. ([57)1 the in- 
tegrand is zero in the bulk phases, therefore the addi- 
tional integral term only renormalizes the interfacial en- 
ergy, and this effect vanishes in the sharp interface limit. 
Nevertheless, the above expression shows that in a stress 
free situation the elastic energy does not vanish in the 
interfacial region, and therefore a mapping of stress to 
temperature is no longer possible as before. 

However, it is intuitively clear that in the sharp inter- 
face limit, where the choice of the interpolation becomes 
irrelevant, we always recover the Clausius-Clapeyron re- 
lation for the considered case of a diagonal eigenstrain, as 
expected. In particular, the stress free branch for the case 
of fixed displacements is indeed a sharp interface result 
and does not depend on the coupling for W ^ (cr//i)^/^. 

VIII. CONCLUDING REMARKS 

In summary, we have developed an analytical approach 
to compute short-range forces between diffuse interfaces 
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in multi-phase-field models based on a mechanical ana- 
log. Even though we have discussed our results in the 
context of grain boundary wetting, the approach is gen- 
eral and should be applicable to a broad class of prob- 
lems. We have found that multi-phase-field formulations 
generally considered in the literature do not describe 
repulsive interactions at large distance when the phase 
fields strictly represent local phase fractions. Motivated 
by this limitation, we have introduced a simple two- 
phase-field model with tunable interaction, which can 
represent both attractive and repulsive interactions. 

This model was only developed here for a bicrystal. 
Therefore, it would be interesting to extend this formu- 
lation to represent an arbitrary number of grains of differ- 
ent crystal orientations, while retaining the flexibility to 
represent both attractive and repulsive interactions be- 
tween different grains. Such a formulation should prove 
valuable to investigate strain localization in the context 
of hot cracking with coupling to elasticity and a shear 
modulus dependent on liquid film width. 

Our study of solute and stress effects has yielded some 
interesting insights that warrant further investigations. 

Firstly, above some threshold concentration, solute ad- 
dition induces coexistence of different grain boundary 
states with different liquid film widths over a finite tem- 
perature range just below melting. This range is gener- 
ally small but is estimated here to increase up to tens of 
degrees for a partition coefficient k ^ 0.01 — 0.1. Interest- 
ingly, the strength of this effect can be understood ana- 
lytically to scale ~ — Infc in the dilute binary alloy limit. 
To our knowledge, this type of bistability has not yet 
been clearly observed in molecular dynamics simulations 
or experiments, although its existence has been predicted 
in other phase-field modeling studies of elemental mate- 
rials and alloys [lllj ESj E^l • light of the present 
results, it would be interesting to test for its existence in 
highly partitioning alloys with interfaces that are struc- 
turally and chemically diffuse. High-energy boundaries 
in such materials appear to be the most likely candidates 
to observe coexistence of equilibrium states with different 
liquid film widths close to melting. 

Secondly, we have found that a uniaxial stress per- 
pendicular to the grain boundary plane is equivalent 
to a temperature change through a standard Clausius- 
Clapeyron relation. Therefore, for repulsive boundaries, 
a large tensile stress of magnitude comparable to a frac- 
tion of the yield stress should suffice to produce an ob- 
servable increase of the liquid film width slightly below 
the melting point. This effect should be readily testable 
by molecular dynamics simulations. 
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APPENDIX A: EQUIVALENCE OF 
CONSTRAINED AND UNCONSTRAINED 
FORMULATIONS OF THE MECHANICAL 
ANALOG 

Here we show explicitly that the elimination of the 
third field leads to the same result as keeping all inde- 
pendent fields and the additional constraint 01 -I- 02 + 
03 = 1. We use a general free-energy density / — 
/(0i, 02, 03,01, 02, 03) and the constraint 0i -t-02 + 03 = 
1. Then the stationary phase field equations are obtained 
from variation of the functional 



F : 



dVf, 



(Al) 



with / — / — Ao(a::) (01 +02 + 03 — 1). Stationarity requires 
SF 



0, i==l,2,3. 



(A2) 



From that we get the expression of the Lagrange multi- 
plier 

and e.g. the first equation 

2^/ _ 9/ _ a/ / 5/ ^ 9/ _ 9/\ 

901 902 903 dx \ 901 902 903 / 



The generalized momenta are 



Pi 



df 



The Hamiltonian is 

- H = F -_pi0, -_P202 -P303- 



(A5) 



(A6) 



Now we can also eliminate the third field from the be- 
ginning, introducing a new free-energy density 

7(01,02,01,02) = 7(01,02,1-01-02, 01, 02, -01-02), 

(A7) 

and we have only two equations of motion 



df d df ^ . ^ ^ 

dpi ax 90J 



(A8) 



They are explicitly 



df 




--( 


^9/ _ 


9/ 


901 


903 


dx \ 


v90i 


903 


9/ 




--( 


^9/ _ 


9/ 


902 


903 


dx ^ 


V902 


903 



0. (AlO) 
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Combining them gives us e.g. the same equation of mo- 
tion (|A4[) above. We can also calculate the energy, and 
define momenta: 



Therefore we get 



Pi 



P2 



df 



a/ _ a/ 

d(j)2 903 ' 



(All) 

(A12) 
(A13) 



and the Hamiltonian 

-H = F-p^4>i-p^^2 = -H, (A14) 
reduces to the same expression as before. Here we used 

03 = -01 - 02- 



and ^ are symmetrical and antisymmetrical about this 
point, respectively. One of the shooting constants can be 
set to an arbitrary value since the problem is invariant 
under a translation along x, i.e. it just fixes the position 
of the origin. Therefore, we integrate from a point far in 
the solid where the asymptotic analytical solutions are 
valid up to the point where d(j)/dx = is reached. We 
then use the remaining shooting parameter to fulfill the 
other boundary condition ■0 = 0. The value of the liquid 
layer thickness can then be extracted by measuring twice 
the distance between the point where ip — —1/2 and the 
endpoint of integration, since the profile is symmetric 
with respect to the latter point. 

For the multi-order parameter model, a similar strat- 
egy can be employed. As definition of the liquid layer 
thickness we use here the distance between the points 
where the two solid fields cross the value 1/2. 

The stable branches can of course also be found by full 
relaxation according to Eq. , but the above procedure 
is more accurate and eflacient also in these cases. 



APPENDIX B: NUMERICAL METHOD 



Here we briefly explain the numerical shooting method 
used to solve the stationary phase-field equations, which 
have the form of a coupled ordinary differential equations 
(ODEs). For simplicity, we describe this method for the 
two-phase- field model of section |V] but the same method 
is applicable to the the multi-order parameter models 
after elimination of one field using the constraint 0i + 

02 + 03 = 1- 

We integrate the ODEs starting in the left grain, i.e. 
for a large negative x with the origin chosen midway be- 
tween the grains. To find out the initial conditions for 
this integration, we linearize the phase field equations 
around the fixpoints, i.e. ip = —l + Sip and — —Scj). 
Then the phase field equations ([^5]) become to first order 
(for simplicity h = a = 1) 



92 



dx-^ 



(8 + 80g)(50 + 6L 



i(8 + 802)<5V 
a 



T-Ti 



M 



Tm 



and only have the exponential solutions 
V' = -l + Ciexp[(8 + 80§)i/2Q,-i/2^- 

= C2exp 



L{T-Tm) \ 
Tm^I ) 



1/2 



which vanish at x — > — cxd, as required for a physically ad- 
missible solution; this requirement fixes one of the two in- 
tegration constants for each second order equation. The 
remaining two constants Ci and C2 are used as adjustable 
parameters in the shooting method to fulfill the bound- 
ary conditions d(j)/dx = and -0 = at the midpoint 
between the two grains, which follow from the fact that 



APPENDIX C: DOUBLE-OBSTACLE 
POTENTIAL 

Instead of the multi-well a multi-obstacle potential is 
often used for phase field simulations 0], and we briefly 
investigate its behavior concerning short-range interface 
interactions here. We refrain from a full analytical and 
numerical treatment and discuss for simplicity only a 
model with a single order parameter. 

The main difference between the double-well and the 
double-obstacle potential is that the latter is defined to 
be infinite outside the the physical regime ^ 7^ 1, and 
this is sketched in Fig. Furthermore, the potential 
has a finite slope at the end points = and = 1. A 
typical choice is 



foo 



h(P{l - 
00 



for 7^ 7^ 1 
else 



(CI) 



instead of the double well potential 

/d^ = 2/102(1 -0)2. (C2) 

Again, we use a standard kinetic term of the type 

fk = (t(V0)2. (C3) 

The central point is now that the double-obstable po- 
tential provides stationary interface solutions which have 
only a finite support, i.e. the phase field differs from the 
trivial values = or = 1 only in a finite region. 
For the double-well potential, the phase field approaches 
these limiting values only exponentially. For the above 
choice of free-energy contributions, the stationary inter- 
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FIG. 27: Mechanical potential energy versus phase-field (f) 
(negative of the free-energy density) for the double-obstacle 
and double-well potentials. 



face solution 4>{x) is given by 



a X- xq < -^7r/2 

+ si\i[{x - Xq)/^]) if - C7r/2 < x ~ xq < ^7r/2 



1 



else 



(C4) 

where xq is the interface position for this onc-dimcnsional 
solution and ^ = {a jK)^!'^ is a measure for the interface 
thickness, as before. 



I Double obstacle 

I Double well 

I 
1 
\ 
\ 








L(T-Tm)^/Ys|Tm 

FIG. 28: Sketch of the liquid layer thickness as function of 
temperature for the double- well and double-obstacle potential 
in the framework of a single-order parameter model. The 
precise location of the curves depends on the model and the 
definition of the melt layer thickness. 



If we repeat the mechanical analog of Fig. [2] for the 
double-obstacle potential, it is clear that the particle tra- 
jectory can only exist for T > Tm since a turning point is 
still present. Therefore the interaction is still attractive. 
The main difference is that the potential does not have 
zero slope near the liquid maximum. The potential has a 
finite slope that does not change as T — Tm 0+. This 
slope implies that the particle has a constant negative ac- 
celeration at the turning point even in this limit. There- 
fore, it spends a finite amount of time near the turning 
point and W does not diverge in this limit. Instead, it 
reaches a maximum value as shown schematically in Fig. 
1281 Exactly at T = Tm, hquid films can exist for any 
W larger than this maximum since the interaction be- 
tween interfaces becomes strictly zero. This also implies 
that the disjoining potential vanishes at a finite W for 
the double-obstacle potential. 
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